R/Importance Sampler.R

In JMacDonaldPhD/Epidemics: Inference For Epidemics (one line, title case)

#' Importance Sampling with known epidemic parameters

#' Epidemic process {x_t ; t \in (0, T)}
#' T is the end of the epidemic
#' x_t is the infectious state of every individual
#' in the population at time t
#'
#' The process is assumed to evolve according to
#' parameters \theta
#'
#' Observed Panel Data {y_t ; t \in (t_1, t_2, ..., t_k)}
#' y_t is the infectious state of a sample of size m from the
#' population

#' y_t are independent conditional on the epidemic process at
#' time t.

# = Distributions =

#' We would like to receive samples from the posterior
#' distribution
#'
#' \pi(x_{0:t}| y_{1:t}) \propto \pi(y_{1:t}| x_{0:t}) \pi(x_{0:t})
#'

#' Do this via importance sampling
#'
#' 1. Propose $x^{i}$ \sim q() (This is a particle)
#' 2. Calculate importance weight \omega(x^{i})
#' (Although this is not possible is cannot calculate \pi(y_{1:t}))
#' 3. Repeat 1-2 N times
#' 4. Calculate Normalised weights \tilde{\omega}(x^{i})
#' (This is possible without \pi(y_{1:t}))
#' 5. Use this weighted sample to estimate integrals involving
#'    posterior distribution.

#' weights \omega are the ratio of the posterior and proposal,
#' measuring the difference telling us how representative the
#' sample is to the posterior. If proposals were made directly
#' from the posterior, all weights would be 1.

#' Diagnostic of the Importance Sampler
#' Effective Sample Size
#'
#' $$ ESS := 1/sum_{i = 1}^{N} \tilde{\omega}(x^{i})^{2} $$
#'
#' How many direct posterior samples the sample
#' obtained through importance sampling is worth
#'
#' If sampled directly from the posterior, all
#' weights would be 1
#'
#'

#' The epidemic process is easily simulated using the Gillespie
#' algorithm. This means using a proposal distribution of the
#' form q() = \pi(x_(0:t)|\theta) is trival. Although this is
#' not really informed by the data
#' = (Can we bias simulation of particle x^{i} according to y?) =

epidemicImportanceSampler <- function(panelData, obsTimes, theta, Nparticles){

  #' 1. Simulate N particles $x^{i}$ \sim q() (This is a particle)


  m = length(panelData[[1]][,2])

  particles = lapply(X = 1:N_p, function(X){
    SIR_Gillespie(N, initial_state = c(rep(1, N - a), rep(2, a)), beta = theta[1], gamma = theta[2], output = "panel",
                  obs_times = obs_times, obs_end = tail(obs_times, n = 1))$panel_data
  })



  #' 2. Calculate Likelihood (which is the main component of the weights)

  if(homogen){
    #' extract panel data from particles
    particle_trans = lapply(particles, transition_data)
    if(length(obs_times) > 2){
      y_given_x = sapply(X = particle_trans, function(X) prod(mapply(extraDistr::dmvhyper, y, X, MoreArgs = list(k = m))))
     }else{
       y_given_x = sapply(X = particle_trans, function(X) extraDistr::dmvhyper(y[[1]], X[[1]], k = m))
     }
  }

  #' 3. Calculate Normalised weights \tilde{\omega}(x^{i})
  ISweights = y_given_x/sum(y_given_x)

  ESS = 1/sum(ISweights^2)


  return(list(ESS = ESS, ISweights = ISweights, particles = particles, y = y))


  #' 4. Use this weighted sample to estimate integrals involving
  #'    posterior distribution.


}